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We discuss how the integrators used for the Hybrid Monte Carlo (HMC) algorithm not only 
approximately conserve some Hamiltonian H but exactly conserve a nearby shadow Hamiltonian 
H, and how the difference AH = H — H may be expressed as an expansion in Poisson brackets. 
By measuring average values of these Poisson brackets over the equilibrium distribution °^ e^^ 
generated by HMC we can find the optimal integrator parameters from a single simulation. We 
show that a good way of doing this in practice is to minimize the variance of AH rather than its 
magnitude, as has been previously suggested. Some details of how to compute Poisson brackets 
for gauge and fermion fields, and for nested and force gradient integrators are also presented. 
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1. Introduction and motivation 

Hybrid Monte Carlo ^ is the algorithm of choice to generate dynamical configurations for 
lattice QCD. This algorithm relies on the introduction of a fictitious momentum for each dynam- 
ical degree of freedom, resulting on a Markov chain with a fixed point exp{—H{q,p)) where the 
Hamiltonian is H = ^p^ + S{q) = T{p)+S{q); ignoring momenta p, we get the desired distribution 

exp{-S{q)). 

The HMC Markov chain alternates two Markov steps: Molecular Dynamics Monte Carlo, 
which consists of a reversible volume-preserving approximate Molecular Dynamics trajectory of 
t/5t steps followed by a Metropolis accept/reject test with acceptance probability mm{l,e^^^); 
and Momentum refreshment from a Gaussian heatbath P{p) °^ e^^ 

1.1 Symplectic integrators 

Symmetric symplectic integrators form a large class of reversible and volume-preserving in- 
tegrators. The idea of a symplectic integrator is to write the evolution operator as exp(T^) = 

exp (^T + ^^1^ = where the Hamiltonian vector field 

dH d dH d , ^ d r^, , . d - ^ 

dq dp dp dq dp dq 

We now make use of the Baker-Campbell-Hausdorff (BCH) formula, which tells us that the 
product of exponentials in any associative algebra can be written as \a{e^l^e^e^/^) — {A-\-B) = 

^ I [A, [A,B]] — 2[5, [A,B]] I H where all the terms on the right hand side are constructed out of 

commutators of A and B with known coefficients. We find that for an STS integrator with step size 
5x the evolution operator for a trajectory of length T may be written as 



exp 
exp 



f + S)5T-^A[S, [S, f]\ + 2[f , [5, f]] ) + &{dT') ' 



24 

T ( f + 5 - ^ f [5, [5, f]] + 2[f , [S, f]]) + ff{5t 



1.2 Shadow Hamiltonians and integrator tuning 

For every symplectic integrator there is a shadow Hamiltonian H that is exactly conserved; 

this may be obtained by replacing the commutators [5, T] in the BCH expansion with the Pois- 
dSdT dSdT ^ 

son bracket {S,T\ = - — r -— — — [EP. For example, the integrator above exactly conserves 

dp dq dq dp 

the shadow Hamiltonian ^57-5 = r + 5 - ^ (^{S, {S,T}} + 2{T, {S, T}}^ 5r^ + ^(St^). We now 
make the simple observation that all symplectic integrators are constructed from the same Pois- 
son brackets (which are extensive quantities). We therefore propose to measure the average values 
of the Poisson brackets ({S, {S, T}}) and {{T,{S,T}}) over a few equilibrated trajectories at the 
parameters of interest and then optimize the integrator (by adjusting the step sizes, order of the 
integration scheme, integrator parameters, number of pseudofermion fields, etc. [Q, ^, |5|] offline) so 
as to minimize the cost. 



2 



Tuning HMC using Poisson brackets 



As a simple example consider the STSTS integrator 
U^TSTsiSzy^'^' = (^e«'^5Tgjf5T^{i-2a)s5Tgif5TgaS5Ty/'" ^j^^^g shadow Hamiltonian is 

6a^-6a + l l-6a,_ r„^ii\o 2 



//srsr5 = ^sr5T5 + {S, {S,T}} + -^{T, {S, T}} j + ^(5t^). (1.1) 

Here we cannot completely eliminate the coefficient of the 0{5x^) contribution as we only 
have one free parameter a, but we can attempt to minimise the cost by adjusting a given the mean 
values {{S, {S,T}}) and ({T, {S, T}}). Naively we could try to minimize the coefficient of in 



(1.1), but we will see below that this is not the best thing to do. 



1.3 Force gradient integrators 

Let us consider again the STSTS integrator, where we set a = g so that the {T, {S, T}} con- 
tribution is eliminated. The remaining leading order Poisson bracket {S, {S,T}} depends only on 
q, which means that we can evaluate the integrator step gi'^ l'^'^}}^^^ explicitly, 

Ufc[5T:)=e(>^e'- e '2 e(>^. 

The force for this integrator step involves second derivatives of the action, and therefore they 
are called Hessian or force gradient integrators [0, |8]]. By putting such an integration step into a 
multistep integrator we can eliminate all the leading G(bz^) terms in A//, as we can see from the 
corresponding shadow Hamiltonian: 

Hfg = T + S- ( 41 {5, {5, {5, {5, r}}}} + 36 {{5, T}, {S, {S, T]]] 

+11 {{SJ},{T, {S, T}}} + 84 {r, {5, {S, {5, r}}}} 
+126 {r, {r, {5, {s,T}}}} + 54 {r, {r, {r, {5, r}}}}) . 



Note that the coefficients of the leading order correction in the shadow Hamiltonian are ap- 
proximately two orders of magnitude smaller than the corresponding coefficients in the Cam- 
postrini integrator [§, 10 1. 



1.4 Nested integrators 

If it is much cheaper to evaluate the force for one part of the action, such as the pure gauge 
part, we can use a nested integrator with a small step size for the inner cheap part. One might 
expect that one could then tune the outer part without reference to the cheap part, but this is not the 
case. 

Let the Hamiltonian he, H = + S\ + S2 with \\S2\\ ^ H^i || and consider a nested integrator 
with a composite step of the form U{8z) = exp ^ (exp ^ exp — exp ^) exp ^ . For the 



inner integrator the BCH formula tell us that ^exp exp exp may be written as 

exp('(5i+f)5T+(a[5i,[5i,f]]+^[f,[5i,f]])^ + ^(5T5) 



3 



Tuning HMC using Poisson brackets 



with oc = and P = j2- Applying the BCH formula again leads to the shadow Hamiltonian 

H = H + (a{S2,{S2,f}} + l5{Sy,{S2,f}} + l5{f,{S2,f}} 

+^(a{5i,{5i,f}} + j8{f,{5i,f}}))5T2 + ^(5T4). 

Observe that the Poisson bracket {Si , {S2,f}} depends on the cheap action but is not supressed 
by any inverse power of m; it is therefore still necessary to measure this quantity in order to optimize 
the integrator. 

2. Computing Poisson brackets 

2.1 Gauge fields 

We must construct the Poisson brackets for gauge fields, where the field variables are con- 
strained to live on a group manifold. To do this we need to use some differential geometry [^]. In 
order to construct a Hamiltonian system on such manifold we need not only a Hamiltonian func- 
tion but also a fundamental closed 2-form CO. On a Lie group manifold this is most easily found 
using the globally defined Maurer-Cartan forms Q' that are dual to the generators and satisfy the 
relation dQ' = — 4cy^0^ A 0*^, where c^^ are the structure constants of the group. We choose to 
define w = -dZi 6'p' = HiiO' A dp' - p'dd') = ^£1(6' A dp' + ^p'c'ji^dJ A 0^): using this funda- 
mental 2-form we can define a Hamiltonian vector field A corresponding to any 0-form A through 
the relation dA{x) = (o{A,x) for all vector fields x. 

For a Hamiltonian of the form H = S + T we find that the leading Poisson brackets that ap- 
pear in the shadow Hamiltonian for a symmetric symplectic integrator are {5, {S,T}} = ei{S)ei{S) 
and {T,{S,T}} = —p'p^eiej{S) where the are the momentum coordinates and the are linear 
differential operators satisfying e, (?7) = TiU for gauge fields U G SU(?i) with generators Tj. 

2.2 Fermions 

Consider the Wilson pseudofermionic action S = (p^^^^ip, and recall that the ei are linear 
differential operators, thus = —^^^^^ei{^)^^^^, and 

p'pjeiej{S) = /jy>^^"' [2ei{^)^-^ej{^)-eiej{^)] ^-^<p. 

ei{^) is straightforward to evaluate given the linearity of the Wilson-Dirac operator in the 
gauge field: we just use Leibniz rule and then replace the gauge field U by TiU . 

3. Results 

3.1 Shadow Hamiltonian and Poisson brackets 

The blue curve in the first plot of figure |I| shows how logjQ \ 5H\ = logjQ \Hf — //,| behaves as 
a function of MD time, compared with the red curve logjg \5H\ for the shadow Hamiltonian up 
to leading non-trivial order in 5t. The simulation uses the STSTS integrator with a = 0.24 and 
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Lattice volume 



(a) Time history of Shadow Hamiltonian. (b) Volume scaling of Poisson brackets. 

Figure 1: Shadow Hamiltonian and Poisson brackets. 



5z = 0.1 for Wilson gauge and fermion actions. This demonstrates that the shadow Hamiltonian is 
indeed conserved. 

The second graph on figure |l| shows how several different Poisson brackets and their fluctua- 
tions depend on the lattice size. As expected the Poisson brackets are more-or-less extensive (they 
grow as L"^); the statistical fluctutations in the Poisson brackets are also shown, and they fall as L^^ 
relative to the mean values as expected. 

3.2 How to tune an integrator? 

We are concerned in minimizing the cost of 
HMC; in our case, this corresponds to maximiz- 
ing the step size 5x while maintaining a reason- 
able acceptance rate. The first step to this goal 
is to find the integrator parameters that maximize 
the acceptance rate for a given value of 5t. Here 
we are going to discuss results for the STSTS inte- 
grator described above, trying to find the optimal 
value for a. 

Omelyan et al. proposed that one should 
minimize (AH^) = (^{H — //)^), as this makes H 
as close to H as possible. However, the amount by 
which AH varies over the equilibrium distribution 
oc e^^ turns out to be considerably smaller than 
the values of AH itself. Therefore, it seems more reasonable to minimize Var(A//), the variance of 
AH over this equilibrium distribution. 

Indeed, figure ^ verifies that (AH) » y^'Vai{AH). If we assume that Hf and are selected 
independently from their equilibrium distributions, which is a goal of HMC, {AH) » (SH) as 
figure ^ also verifies. We can also conclude that the initial and final distributions seem to be 
equivalent — of course, Hf is not distributed according to the equilibrium distribution as //, is, but 
its distribution does not differ significantly. 



Figure 2: Histogram of AH at the start (blue) 
and end (red) of the trajectories. 
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-logio |<5H>| 

-loglO <6HA2 / 2> 
-loglO <1 - min(l,eA(-5H))> 
-loglO Var(AH) 
-loglQ <AHA2/1QA4> 



(a) STSTS integrator. 




(b) Two level STSTS. 



Figure 3: Tuning plots. 



Function 



\{5H)\ 0.16749 

{5H^/2) 0.17765 

Var(A//) 0.18260 

(l-min(l,f'-^^)) 0.18664 



0.24952 



In figure [3(a]| , we see plots of several quan- 
tities, besides Var(A//), we could minimize to 
optimize the integrator. The curves were com- 
puted using the Poisson brackets computed at a = 
0.24, whereas the red points are measurements of 
(5//^/2) at different a values. The good agree- 
ment between the measured and predicted loca- 
tion of the minimum gives us confidence that we 
can find the correct behaviour of the quantities 
of interest by measuring the Poisson brackets at 
a single value of the integrator parameters. 

In table | we can see the optimal a values for the quantities considered. We see that the 
minima for \{dH)\, {dH^/2) and (1 — min(l,e^^^)) are close to the minimum of Var(A//). 

Figure |3(b)| shows similar results for tuning the parameters for a dynamical fermion computa- 
tion on a 8"* lattice with a Wilson gauge action with j8 = 5.6 and Wilson fermions with k: = 0.1575. 
Here we minimize {5H^). We used a two level STSTS integrator with two gauge steps per fermion 
step, and a trajectory length of one. The yellow point shows values of the a parameters at which 
the Poisson brackets were measured. 



Table 1: Optimal values for a. 



3.3 Force gradient integrators 



In this subsection, we show results for the force gradient integrator defined in section 1.3 
obtained with the Wilson gauge action at j8 = 5.6 on a 4^* volume, comparing with a second order 
Omelyan integrator (figure Note that the scaling for the force gradient integrator (black data in 
figure |4(b)| ) is much better than for the Omelyan integrator (green data in figure 4{b] ). 
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Comparison of 2nd order Omelyan and force Gradient STSTS ^^.^ij^g (he force gradient integrator 

Wilson gauge action, V-4 , G=5.6, 5t=0.05 ,. .-t ,, ^ , 

° " ' • ■ Wilson gauge action. , p=5.& 

I ' \ ' \ ' \ ' \ ' I b ^ ^ ^ 

I _ ; 




Trajectory Step size 

(a) MD history. (b) Scaling. 



Figure 4: Results for the STSTS force gradient integrator (green data), compared with data for second 
order Omelyan integrator (black data). Note also the data for the shadow Hamiltonian of Omelyan integrator 
(red data). 

4. Conclusions 

We have shown that a good strategy to optimize HMC integrators is to minimize the variance 
of AH over the equilibrium distribution e^^ , rather than minimizing \AH\ itself, as was previously 
proposed. We have outlined how the Poisson brackets required to compute AH may be evaluated for 
gauge theories and systems with dynamical fermions. We have also carried out initial investigations 
with nested integrators and force gradient integrators. We hope to present more details of our 
results, and data for more realistic computations soon. 
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